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We investigate the use of hybrid techniques in complex processes of infectious diseases. Since pre- 
dictive disease models in biomedicine require a multiscale approach for understanding the molecule- 
cell-tissue-organ-body interactions, heterogeneous methodologies are often employed for describing 
the different biological scales. Hybrid models provide effective means for complex disease mod- 
elling where the action and dosage of a drug or a therapy could be meaningfully investigated: the 
infection dynamics can be classically described in a continuous fashion, while the scheduling of mul- 
tiple treatment discretely. We define an algebraic language for specifying general disease processes 
and multiple treatments, from which a semantics in terms of hybrid dynamical system can be de- 
rived. Then, the application of control-theoretic tools is proposed in order to compute the optimal 
scheduling of multiple therapies. The potentialities of our approach are shown in the case study of the 
SIR epidemic model and we discuss its applicability on osteomyelitis, a bacterial infection affecting 
the bone remodelling system in a specific and multiscale manner. We report that formal languages 
are helpful in giving a general homogeneous formulation for the different scales involved in a mul- 
tiscale disease process; and that the combination of hybrid modelling and control theory provides 
solid grounds for computational medicine. 



1 How many scales does it take to describe an epidemics? 

In this work we investigate multi-methodology approaches for the modelling of complex infectious dis- 
eases and for finding the best strategy of intervention, able to avoid or limit the spread of the disease 
and its effect on the affected organisms. Infectious diseases can be broadly classified into three groups 
(Fig. [TJ: acute, latent persistent and chronic persistent. Acute diseases like the common cold, the Rhi- 
novirus, the Yellow Fever, the Influenza or some strains of the Staphylococcus Aureus (Osteomyelitis) 
are characterized by a single disease episode after which they do not occur anymore. Latent persistent 
ones like the Herpes simplex, the Varicella-zoster or the Measles-SSPE arises also after the first disease 
episode and a non-infectious and latent period. Chronic persistent diseases (e.g. Hepatitis B, HIV, HTLV- 
1 leukemia, chronic Osteomyelitis) can protract their effects on the host organism for several years. 

Epidemic modelling is one of the most established tools for predicting the progress and the spread 
of a disease in large populations, which are commonly divided into compartments: Susceptible, Exposed, 
Infected and Recovered. According to the possible flows among such compartments, different variants 
arises, like the SIR (Susceptible^Infected^Recovered) model HU, SIS, SIRS, SEIR and SEIRS. Fur- 
ther compartments can be taken into account like the population of immune infants, usually indicated 
with M. Similar compartments can be identified not only at the human population level, but also at the 
cellular population level, where the pathogen typically acts by infecting susceptible cells. As a matter of 
fact, infections are characterized by multiscale dynamics that affect the organism at multiple levels in the 
biological hierarchy, as shown on the left side of Fig. [2] at the intracellular level, in case of infection of a 
cell by a pathogen; at the intercellular/cellular population level, in case of infection of susceptible cells 
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Common cold, Rhinovirus, Infectious cancers, Herpes simplex, Varicella-zoster, Hepatitis B, HIV, S. aureus (Osteomyelitis), 

Yellow Fever, Influenza, S. aureus (Osteomyelitis) Measles-SSPE HTLV-1 leukemia, Prion (creutzfeldt jakob) 

Figure 1 : Simplified classification of infectious diseases. The disease intensity follows a Gaussian curve 
approximately partitioned in three stages. In the first one, the infection starts but any symptom is re- 
ported. In the second one called "disease episode" symptoms are reported and the disease can be trans- 
mitted to other individuals. In the last stage, the disease is over if it is acute or latent persistent. However 
some diseases can become chronic. 



by infected cells; at the tissue level, when multi-cellular ensembles are involved; at the organ/individual 
level, when infection spreads to other parts of the body; and at the human population level, if the disease 
is transmitted among individuals. A further scale above the human population one could be the con- 
sciousness level for which several psychological and cognitive behaviour models exists. Just consider 
that psychological states like fear or stress can even affect (negatively) the immune system: a demonstra- 
tion that all the scales in the biological hierarchy are intimately connected to each other. Their inherent 
multiscale dynamics make the modelling of infective processes a challenging and intriguing field of re- 
search. Indeed, pathogens and infective agents generally act at the cellular scale, but they need to exploit 
mechanisms at the human population level to spread and survive in other host organisms. 

1.1 Formal languages, hybrid modelling and control in multiscale infectious diseases 

We believe that hybrid modelling could help in unravelling the complexity of the multiscale dynamics 
occurring in such diseases. We refer to 'hybrid modelling' not just with its classical meaning, i.e. the 
modelling of those systems characterized by the co-existence of continuous and discrete dynamics, but 
also with a methodological meaning. In multiscale systems different scales are typically approached 
with different methodologies, thus leading to the co-existence of heterogeneous modelling techniques, in 
other words to a methodologically hybrid modelling approach. The cellular level is typically represented 
with agent-based or ordinary differential equation models, while the tissue and organ levels are often 
described using image-based finite element modelling (partial differential equations). In the context 
of disease modelling the determination of suitable intervention strategies, including drug and therapy 
administration, adds another level of description that affects multiple biological scales. 

This kind of integrative models are composed of single-scale models, describing the biological pro- 
cess at different characteristic space-time scales, and scale bridging models, which define how the single- 
scale models are coupled to each other ll22ll . Higher level languages and formalisms can help in giving a 
general homogeneous formulation for the different scales in a multiscale biological system. In general, 
that formal description cannot be directly executed or simulated, but it is able to support multiple seman- 
tics (e.g. transition systems, differential equations and Markov chains). In particular process-algebraic 
languages are a formal notation initially developed for modelling software systems, but in the last decade 
it has been extensively used and extended in order to describe biological systems. In a seminal paper 0, 
Cardelli showed that a subset of stochastic CCS is powerful enough to encode both systems of reactions 
in the stoichiometric form and systems of ordinary differential equations. Additionally, a recent work of 
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Figure 2: Modelling of multiscale infective dynamics. Infections involve the intracellular level (a cell 
infected by a pathogen); the intercellular level (infection among cells); the tissue level (infection among 
collections of cells); the organ level (infection among different parts of the body); and the human popu- 
lation level (transmission among individuals). Each biological scale can be implemented with different 
modelling techniques like Gillespie's SSA, Continuous Time Markov Chains (CTMC), ordinary/partial 
differential equations (ODE/PDE), agent-based models (ABM), or finite state machines (FSM). Such 
heterogeneous semantics can be wrapped by a common formal language. 



the authors |[T8l[T9l demonstrates the effectiveness of a hybrid approach where a process-algebraic spec- 
ification level is translated into a runnable stochastic agent-based model in the bone remodelling case 
study, while in [2 ] the same case study is approached with different semantics according to the biological 
property to analyse. 

Consequently, formal languages could represent a wrapping language able to homogeneously de- 
scribe the different scales of a complex multiscale system, where each level can be instantiated into 
a runnable model according to the most suitable semantics. Figure [2] sketches the idea of the hetero- 
geneous semantics associated to a multiscale model of infection through a common formal language. 
Furthermore, the importance of being hybrid in the classical sense is demonstrated by several biologi- 
cal evidences. For instance, genetic regulatory networks naturally exhibits hybrid dynamics, for which 
the continuous concentration of proteins is interrupted by the discrete switches dictated by changes in 
gene expression 12 .]___]• Moreover in population biology it happens that some species are present in 
high concentration, so they can be approximately modelled by continuous variables; conversely, small 
populations are opportunely modelled as discrete stochastic variables |[T4l . 

Another crucial aspect is related to the self -regulation and control of biological systems: in normal 
conditions, biological entities and functions are self-regulated and several multiscale control mechanisms 
naturally exists. One of the most striking example is the control operated by the immune system that 
protect the organism in case of disease. In fact, lymphocytes are able to detect the presence of a pathogen 
and produce an appropriate immune response by secreting immunoglobulins. However, in many cases 
of severe diseases the immune system cannot apply an effective control anymore. Drugs and therapies 
represent a form of control which is external to the organism and is extremely significant when dealing 
with models of diseases that can be limited by appropriate medical interventions. The control operated by 
medical therapies can be much more effective but needs to formulate control laws that take into account 
their (negative) impact on the organism due to possible side effects. While the drug administration 
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strategy is determined at the human population level, the therapy affects multiple lower scales (in most 
of the cases, the cellular and subcellular scale). Other examples of multiscale control are the mechanical 
control at the organ and tissue level that regulates the functioning of bone cells in the bone remodelling 
process; or in an infection scenario, the cognitive behaviour at the human level that avoids or limits the 
interactions with other individuals possibly affected by the disease. A work by Bagnoli et al. HI embeds 
this aspect of human behaviour into a SIR epidemic model by allowing individuals to perceive the risk 
of being infected by ill neighbours. Therefore, the use of control-theoretic tools seems promising for 
describing the biological mechanisms of multiscale self-regulation and self-adaptiveness, as well as for 
implementing externally imposed control mechanisms. 

In this work we present a general model for complex multiscale infectious diseases and for scheduling 
optimal medical treatments. We employ a hybrid modelling approach both in the classical and in the 
methodological sense. Indeed, we formulate the disease process by means of a high-level process algebra 
called D-CGF, from which a hybrid (in the classical sense) semantics can be derived. This semantics is 
given in terms of a hybrid dynamical system, where the cell population and the infection dynamics are 
described in a continuous fashion, while the dosage of multiple therapies is implemented as discrete 
switches. Then, we propose the application of model predictive control (MPC) tools for computing the 
optimal scheduling of multiple therapies. The potentialities of our approach are shown in the case study 
of the SIR epidemic model and we discuss its applicability on the case study of osteomyelitis, a bacterial 
infection affecting the bone remodelling system. 

The paper is organized as follows. Section [2] presents the algebraic language D-CGF for generic 
disease processes and its semantics in terms of hybrid dynamical systems, by means of the SIR example. 
Section [3]reports the optimal controlled solutions for scheduling multiple therapies in the SIR model and 
show how this approach can be useful to the osteomyelitis case study. Discussions and conclusions are 
given in Section [4] 

2 General process-algebraic formulation of infectious diseases 

In this part we propose a generic formulation of complex infectious diseases, by defining a variant of 
the Chemical Ground Form (CGF) stochastic process algebra. Stochastic process algebras extend clas- 
sical ones with quantitative information in form of action rates and have associated a Continuous Time 
Markov Chains (CTMC) semantics from which a continuous and less computationally demanding ap- 
proximation in terms of ODE can be derived. The multiplicity of semantics it supports has made formal 
process-algebraic languages one of the most used tools in biological modelling. A work by Cardelli Q 
shows how several encodings can be implemented among specifications in the CGF algebra (a subset of 
stochastic CCS), systems of reactions in the stoichiometric form and ODE systems. Bortolussi and Poli- 
criti (H extend this work by defining a semantics in terms of hybrid automata [11], defining the requisites 
according to which a set of processes can be treated as discrete control states, while others have the usual 
continuous interpretation. Several other methods for associating hybrid semantics to quantitative process 
algebras have been proposed, among which |9j |3 . 

The formal modelling language employed in this paper is called Disease Chemical Ground Form (D- 
CGF), and it is a variant of CGF for describing complex disease processes. Differently from fPTll . where 
a stochastic process algebra and its continuous semantics are used to describe a classical epidemic model, 
we attempt to provide a framework for more generic diseases and, although the syntax of D-CGF does 
not radically differ from that of CGF, we define a novel constructive procedure for deriving a semantics 
in terms of hybrid dynamical systems, by distinguishing the following two kinds of processes: 
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Table 1: Syntax of D-CGF. 



• Individuals: standard CGF processes that belong to some species and that collectively represent the 
population in the disease model. Different species can describe for instance the various compart- 
ments of an epidemic model (susceptible, infected, and etc.), a pathogen, or a particular mutation. 
Individuals are interpreted as continuous variables in the hybrid semantics. 

• Therapies: processes for modelling interventions on the disease scenario. They could represent 
for instance the dosage of a drug or a change in the environment that discretely alter the disease 
dynamics. Based on the hybrid semantics of CGF |6], therapy processes are subject to some 
restrictions in order to be interpreted as the discrete switches of the hybrid dynamical system. 

D-CFG will be illustrated by means of an epidemic example, even if this language is suitable to describe 
also models of virus infection, infectious cancers or other kinds of diseases where discrete intervention 
policies has to be modelled. 

2.1 The Disease Chemical Ground Form 

Here we present the Disease Chemical Ground Form (D-CGF), a variant of the Chemical Ground Form 
(CGF) [8 ]) targeted to describe complex models of disease in a process-algebraic fashion. The syntax of 
D-CGF is given by the grammar in Table [T] A D-CGF model is given by a set of species S, an initial 
population P, a set of therapies T, and an initial combination of therapies C. The first notable difference 
from the algebra CGF is that in D-CGF models two disjoint sets of processes are distinguished: 

• the set of species (a set of individual species definitions), the individual species and the population 
(the parallel composition of individuals); and 

• the set of therapy definitions, the therapy and the combination of therapies (the parallel composi- 
tion of therapies). 

Individual species and therapies are defined by the alternative composition (+) of action-prefixed terms. 
Actions 71 G IT are indexed by a quantity r corresponding to the rate of the exponential distribution that 
determines their duration, and are the usual ones: x r (internal action), ?x r (input action) and \x r (output 
action). The motivation for separating these two sets of processes is related to the hybrid semantics 
of D-CGF: species are interpreted as continuous variables, while therapies as discrete control states. 
In addition we will need to put some restrictions on the therapy processes, for formally justifying this 
separation. In the remainder of this section, a simple epidemic model will serve as the running example 
for illustrating our approach. 
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S^ h s + s 

I — >j I + S Birth of susceptible 
R^ b R + S 



S — >■» Death of susceptible 
I — > M Death of infected 
R^u® Death of recovered 



S + 1 — » ^ I + 1 Infection of susceptible 



/ — ^ v R Infected becomes recovered 



Table 2: Reactions in the SIR epidemic model with open population and no disease-dependent death rate. 
2.1.1 D-CGF by example 

In this part, we illustrate some of the technical features of the algebra by means of an epidemic example: 
the SIR (Susceptible— > Infected— > Recovered) model lfl2ll . It describes the dynamics at the epidemiolog- 
ical level of a population consisting of individuals susceptible to the disease (S), those infected (I) and 
those recovered (/?). Here we consider an open-population (i.e. with births and deaths) SIR model, given 
by the list of reactions in Table [2j 

According to the conversion rules for CGF, we can translate this list of reactions in a D-CGF model 
^ = (5,P,0,0) with no therapies. 



where Pq is the initial population. Given an action % G IT, we denote with react (n) the multiset of species 
consumed by %, with prod{%) the multiset of those produced by 7T. #(X,P) denotes the number of X 
occurring in the population P. A(n,X) = #(X,prod(n)) — #(X, react (n)) denotes the net variation of 
X due to %. For instance in the above example, reactii) = {S,I}, prod{i) = {1,1}, A(i,S) = —1 and 



The procedure for extracting a system of ODEs requires to build the so-called stoichiometric matrix 
M, which has one row for each species and one column for each action. A cell M[X,7i] takes the net 
variation of X due to %: A(tt,X). Then a rate vector </> has to be defined in the following way. Let r{%) 
denote the rate of action %. Then for each action % € II, 



Finally, the differential equations are given by X = M ■ <p . The stoichiometric matrix, the rate vector and 
the resulting ODEs for the SIR example are shown below. 



S= Tl r {S\\S) + t£.0 + ift.I 
R= t*.(*||S) + t£ 2 .0 



P : Po, 



A(i,7) = 1. 




if react(%) = 
if react(n) = {X} 
if react (n) = {X,Y} 
if react{n) = {X,X} 
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S + Tl on — ¥ p R + Tl on Susceptible becomes recovered with therapy 1 

I + T2 „ — R + T2 on Infected becomes recovered with therapy 2 

T^off ~^ri m T\ on Therapy 1 switched on 

Tlon ~~ M // Tloff Therapy 1 switched off 

T2 ff — > r 2 B „ T2 un Therapy 2 switched on 

T2 on —>r2 a ft T2 ff Therapy 2 switched off 

Table 3: Additional therapy-specific reactions in the SIR epidemic model with open population and 
medical treatments at Table |2l 



Now we modify the SIR example in order to include two therapies Tl and T2 that makes immune 
the susceptible population and increase the mutation rate from infected to recovered, respectively. Tl on 
(T2 0n ) and T\ ff (T2 on ) denote the therapy being administered or not, respectively. Therefore, Tl can 
be thought as a vaccination that makes immune the susceptible individuals, while T2 as a therapy for 
combating the course of the disease. The therapy- specific reactions introduced are listed in Table [3] In 
this case, the associated D-CGF model = (S,P,T,C) has been extended with the therapy processes 



and is defined as follows: 



S = t|j .(S\\S) + t£ .0 + ?fi .1 + ? / .R 

I = t|.(7||S) + t£.0 + li^.R + {fid + 7h k .R P : P 

r= 4a r \w + **2-° 

Tl on = If. Tl„, + x r lf f .T\ off 
T2 off= 

T2„„ = \h v+k .T2 m + & f .T2 n 



rioff 

loff '■"ff C ■ T\ ff II T2 f f 



hoff 1 z °ff 



2.2 Hybrid semantics of D-CGF 

In the following we will give the definition of switching therapy (ST), and show that STs constitute the 
discrete switches in the hybrid dynamical system semantics where the different species are the continuous 
variables. The definition of switching therapy is broadly inspired by the notion of Control Automata 10, 
i.e. the control structure that can be identified in the hybrid automata-based semantics of a CGF model. 
We also provide a constructing procedure for determining the STs associated to a D-CGF model and 
therefore its semantics in terms of hybrid dynamical systems. 

In order to interpret the set therapy terms T as a set of discrete switches, we need to identify collec- 
tions of terms in which exactly one term is active in every combination of therapies reachable from the 
initial one. We call such a collection a Switching Therapy (ST). The intuition is that a switching therapy 
models a discrete component whose terms represent its internal states, since exactly one of them must 
be active at each step. We give the definition of ST and of well-formed therapies for imposing additional 
restrictions on the construction of therapy terms in order to ensure a correct discrete interpretation of 
them. 

Definition 1 (Switching Therapy) Let ^# = (S, P, T, C) be a D-CGF model, and 7} C T a set of therapy 
terms. T\ is a Switching Therapy (ST) if the following conditions hold: 

1. Exactly one of the terms U £ 7] is active in the initial combination of therapies: #(7],C) = 
LueT i HU,C) = \. 
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2. Each action must conserve the concentration of terms in 7} and cannot involve more than one 
reagent in 7}: Vft G II. (# (7J-, react(n)) = #{Tj,prod{tt)) < 1). 

J. For eac/z action % G II smc/z ?/j<3? U\ G react{n) and U2 G prod{%), with U\,U2 &Ti, Ui ^ U2 ==>• 
(react{%) = {U{\ A prod(n) = {U 2 }). 

Conditions 1+2 ensure that exactly one term of the ST is active in every reachable configuration, since 
condition 1 requires that the initial concentration of 7]-terms must be equal to one and condition 2 tells 
that such concentration is conserved. Additionally, condition 2 implies that there are no actions able to 
modify the concentration of 7}-terms and thus that the internal state of the ST 7} can be changed only by 
actions having a 7}-term in its reactants. Finally, condition 3 requires that every action causing the switch 
of a 7}-term from U\ to U2 must be an internal action of U\, of the form U\ = ...+ Z.U2 + . . .. 

Definition 2 (Well-formed therapies) Let Jt = (S, P, T, C) be a D-CGF model andletC = U\\\...\\ U n , 
with n > 1. T is a set of well-formed therapies if there exists a partition 37* = {T\,...,T n \ofT such that 
7J G 3~ is a switching therapy and £/, G T, for all i = 1 , . . . , n. 

We present a method for extracting the switching therapies and show how the semantics in terms hybrid 
dynamical systems can be constructed. Given a D-CGF model = {S,P, T,C), the procedure exploits 
the stoichiometric matrix M, with which it is easy to check the conservation of switching terms and 
condition 3 of Definition [TJ Some necessary conditions has to be formulated on M for ensuring the 
well-formedness of T. 

Definition 3 (Necessary conditions for well-formedness) Let ^# = (5, P, T, C) be a D-CGF model and 
M be the associated stoichiometric matrix. Then, the following conditions are necessary for the well- 
formedness of T 

1. The elements of M restricted to T take values in { — 1,0,1}." V7T G IT. \/U G T. M\x[U, 7t] G 

{-1,0,1}. 

2. Conservation of T -terms: Vtt G IT. Lj/ e r M[t/, k] = 0; 

3. Exclusive switch (I): Mil G IT. 3<i£/ G T. M[U,7t] = —1. This condition and the conservation 
condition additionally imply that it exists exactly one V G T such that M[V, 7t} = \. 

4. Exclusive switch (2): Vft G IT. (3U G T. M[U,x] = -1 (VX G S.M[X,n] = AK internal)). 
The stoichiometric matrix M associated to the SIR example with therapies is given below. 
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It is easy to check that the conditions listed above are met in M.Such conditions on the stoichiometric 
matrix are not sufficient to ensure the well-formedness of the set of therapies T. Indeed, Def. [3] corre- 
sponds to stating that the whole set T meets conditions 2+3 in the definition of switching therapy, but it 
is not possible from the only stoichiometric matrix to prove also condition 1 (exactly one term in T is 
active in the initial configuration). Moreover recalling that T is well-formed if it can be partitioned into 
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a set of switching therapies, such necessary conditions are defined on the entire set of therapies T, thus 
considering only the trivial partition {T}. 

The algorithm for extracting the STs consists in building a graph 'S = (T,E) called ST-graph, where 
T is the set of process terms and arcs connect couples {Ui^Uf) of therapy terms that are involved in a 
switch, i.e. M\U\,n\ = -1 andM\ T [U2, It] = 1 for some %. 

Definition 4 (ST-graph) Let ^ = (S,P,T,C) be a D-CGF model and M be the associated stoichio- 
metric matrix. The ST-graph ( Sj( = (T,E) associated to ^# is a directed graph whose vertices are the 
T -terms of Jg and E = {(U U U 2 ) G T x T \ 3n G n. (M[U h x) = -1 AM ]T [U 2 ,7c] = 1)}. 

Now it is possible to check the well-formedness of T by considering the set of connected components 
of ^tySjt) (see Proposition [5]). In particular, if the stoichiometric matrix M meets the necessary 
conditions in Def. [3]and the set of connected components ^(^S^g) = {G\ , G„} is such that for each 
G; there is exactly one term active in the initial configuration, then T is well-formed and the nodes of 
each Gj form a switching therapy. Note that the sets of nodes of the connected components in an ST- 
graph form, by definition, a partition of T . This is necessary to guarantee that, according to Def. [2] a set 
well-formed therapies can be partitioned into a set of switching therapies. In addition, the therapy terms 
in each G, meet the Condition 2 of Def. [T] about the conservation of terms in a switching therapy, because 
the connected components G, are (of course) mutually disconnected, i.e. there cannot exist any arc in the 
ST-graph connecting them that is, no terms can flow between them. 

Proposition5 Let jM = (S,P,T,C) be a D-CGF model, M be the associated stoichiometric matrix, 
< ^ 7 (^-#) = {Gi = (T[,Ei),. . . ,G„ = (T n ,E n )} be the set of connected components of the associated ST- 
graph and ST = {T\ T n }. T is well-formed and each Ti G S is a switching therapy if the 
following conditions hold: 

1. M meets the necessary requirements at Def. [?]■ 

2. for each Tj G 2T, exactly one of the terms in 7} is active in the initial combination of therapies: 
#(7J,C) = 1; and 

3. for each T; £ J? and for each action % G IT, the number ofTi-terms in react{n) is at most one: 
#(7], react (it)) < 1. 

The ST-graph for the SIR example is given in Fig. [3] (a). It is straightforward to check that its 
connected components are such that exactly one term in each G, is active in the initial combination of 
therapies C = T\ Q ff \ T2 a ff, and that for each Gj there is no action % whose reactants involve more 
than one term in G ( . Therefore the set of therapy definitions T of the SIR example is well-formed and 
5? = {{Tl ff,Tl on },{T2 ff,T2 on }} is the set of switching terms. 

The combination of the different switching therapies allows us to extract the possible discrete modes 
in the hybrid semantics. We define a mode graph structure as the Cartesian product of the switching 
terms in a ST-graph. Figure|3](b) shows the mode graph for the SIR example. 

Definition 6 (Mode graph) Let = (S,P, T,C) be a D-CGF model and c g{<$^) = {G b . . . , G,J be 

the set of connected components of the associated ST-graph. The mode graph MG j( is defined by the 
Cartesian product G\ x . . . x G„. 

Note that transitions between modes could have associated exponentially distributed delays determined 
by the stochastic rates of the algebraic specification. In our context, we omit them assuming that such 
transitions are instantaneous, also because discrete modes will represent the control inputs of the hybrid 
dynamical system semantics. 
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(a) (b) 



Figure 3: ST-graph for the SIR example (a) and its corresponding mode graph (b). The connected 
components of the ST-graph are such that there is exactly one term in each G, active in C and for each 
Gi, no action % involves more than one G;-term in react{%). This ensures that T is well formed and 
the connected nodes are switching therapies. In (b) four distinct discrete modes are generated from the 
Cartesian product of G\ and G2. 

As pointed out also in [6], Petri Nets could have been similarly applied in order to study the con- 
servation properties in a switching therapy. In particular, a set of switching therapies could be seen as a 
set of strictly conservative (i.e. constant number of tokens) and 1-bounded (i.e. the maximum number of 
token is 1) Petri Nets. 

2.2.1 Hybrid dynamical system of D-CGF models 

Here we describe the semantics of a D-CGF model in terms of a particular class of hybrid dynamical 
systems, called controlled switched systems (CSS) iTToTl . While in the general formulation of hybrid 
dynamical system the external control input and the discrete operation mode are distinct, in a CSS the 
external controller produces a switching signal (i.e. the discrete mode) that is given in input to the plant 
(i.e. the controlled system). The basic form of a controlled switched system is the following: 

x=f(x,q) 
y =g(x,q), 

where x G W is the continuous state, f is a vector- valued smooth function, g is the output function and q is 
the discrete operation mode. The above equations are often written with the form x = f q (x) y = gq(x), 
stressing the different dynamics and outputs under different operation modes. Note that a controlled 
switched system can also be seen as a hybrid dynamical system with controlled discrete inputs. Figure [4] 
shows the control loop in a CSS. We show how to derive a CSS from a D-CGF model, by following 
the method elaborated by Cardelli for the CGF algebra as regards the continuous part of the semantics, 
but considering also the switching therapies and the mode graph as regards the discrete part. Recalling 
that each node of a mode graph MGjt = ({q\ = {U\\, . . . ,Ui n }, . . . = {t/jti, . • . , U^}}, E) models a 
particular combination of therapies, a mode qi being active means that therapies U G qt are active and 
that those in T \ qi are not. Therefore a modified version of the rate vector will be used: 

$* = <!>[%,%}, Ueq u VeT\q u 

which is obtained by zeroing all the elements of that contain a therapy term not belonging to qt, and 
omitting (substitution with 1) those in q\ . The resulting controlled switched system would be: 

x = f ? ,(x)=M| S -0 ?i , 
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Figure 4: Control loop in a controlled switched system. The plant defines the piecewise smooth dynamics 
of the system x and the observable output y. The output could represent patient's medical data obtained 
after a visit. Depending on y, the controller acts as therapy scheduler and determines the operation mode 
q (e.g. a particular combination of therapies) of the plant. 



where q, G MGji is a node of the mode graph, and Mw * s tne stoichiometric matrix restricted to the set 
of species S. The operation modes q, the modified rate vector </> q and the resulting CSS x of the SIR 
example are shown below. 
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r S = bN-pSI-llS 
Mllj^l l = PSI-fU-Vl 
I R = vI-nR 



( S = bN-pSI-fiS-pS 
I R = vl- y,R + pS 



( S = bN-pSI-llS 
x(?3) = < i = PSI-pI-(V + k)I x(q 4 ) 
{ R=(v + k)I-LiR 

Alternatively, we could consider the equivalent formulation which exploits the sets of switching terms 
{Tl ff, Tl on } and {T2 ff, T2 on } as discrete control variables in the following way. 



S = bN-pSI-tiS-pS 
i = fiSI-iil-{v + k)l 
R = (v + k)I-fiR + pS 



if Tl on e qi 
ifTl off eqi 



T 2 (qi) = 



if T2„„ 6 qi 
if T2 off e H 



S = bN-PSI-llS- T\ {qi)pS 

t = pSI-nI-(v + T 2 ( qi )k)I 

R = (V + T 2 ( qi )k)I -iiR + n ( qi )pS 



Note that this form is applicable because in the SIR example each switching term contains exactly two 
terms and can consequently have assigned a value in {0, l}. 



3 Hybrid semantics and optimal control of infectious processes 

In this section we define the optimal control law for therapy scheduling in the SIR model by means of 

Model Predictive Control (MPC) [4], and we report the optimally controlled solutions with parameters 

that are typical of the measles |[20l . A similar problem has been addressed in 11211 for the scheduling of 

hormone therapy in a model for prostate cancer. 

In order to cope with the nonlinear dynamics of the SIR model, the optimal scheduling of multiple 
therapies is computed by solving an on-line MPC problem in the Multi-Parametric Toolbox (MPT) lfl3l . 
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The optimal control law can be formulated by means of the following constrained finite-time optimal 
control (CFTOC) problem, where the original continuous-time dynamics is simply discretized with the 
Euler method: 



l + (b-pi(k)-n)At 





bAt 

1 + (pS(k) - /i - v)Af 
vAt 



bAt 


1 -flAt 





= -pS(k)At 







"l 





" 


B = 





-kI(k)At 


Q = 





10 







pS(k)At 


kI(k)At 










0.5 



™"£j=olltfu(*)lli + H0x(fc)||i 

subj. to x(k) e [0, l] 3 
x(T) £ X se t 

x{k + At) =Ax(k)+Bu(k) 
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where At = 1 /365 is the discrete time step corresponding to one day; u(k) = [Tl(k) T2(k)] T € {0, l} 2 
is the discrete input (the combination of the two therapies) and x(k) = [S(k) I{k) R(k)] T € [0, l] 3 is the 
systems' continuous state vector at time k; T = 3 is the prediction horizon (corresponding to three days); 
SC set is a polytope identifying the set of terminal states; ||x||i = |5| + |/| + \R\ and ||u||i = |Z1| + \T2\; 
R and Q are the weights on the controlled inputs and on the system state, resp. The terminal set 3C set 
represents an invariant on the final states and can be described by the convex hull of its vertices: X set = 
Conv({ [1 0] T , [0 l] r }). In this way, the terminal set contains all the states [S I R] T such that / = (no 
infected individuals), S,R > and S + R = 1 (we do not exceed the population limit). Matrix Q has been 
set so that the state variable Infected has associated a weight equals to 10 and thus is slightly penalized 
over variables Susceptible (weight 1) and Recovered (weight 0.5). Finally, we assume for simplicity that 
the internal state can be observed and we set the output vector y to the state vector x. 

Matrix R allows us to implement different therapy scheduling strategies. The higher R, the higher 
penalty is given to drug dosage. In particular if R is assigned a large value, then the strategy is to avoid 
drug dosage as much as possible. This can be the case of a therapy with severe side effects or of a 
patient in an early stage of the disease that can be treated even with a small dosage. On the contrary, a 
low value to R leads to a strategy where drug dosage is much more prominent (e.g. little side effects or 
mature stage of disease). Additionally we can assign different weights to different drugs according to 
their therapeutic impact on the patient, thus possibly enabling a wider spectrum of control strategies. We 
computed the optimally controlled solutions considering three scenarios obtained by varying the weights 
on the controlled therapies: 



1. Low penalties to Tl and T2: R 



2. High penalty to Tl: R 



3. High penalty to T2: R 







100 
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100 



0.1 
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Results reported in Figure [5] tell that solutions under scenario 1 and scenario 2 are the same (Fig. [5] 
(a)), suggesting that therapy T 1 does not need to be administered in the optimal strategy of treatment in 
order to fulfil the terminal constraints. Indeed Fig.[5](b) shows that when T2 has a high weight (scenario 
3), no control moves are performed, so indicating that Tl is ineffective with respect to T2. 



3.1 Osteomyelitis, a truly multiscale infection 

Osteomyelitis is a a bone pathology caused by bacteria infection (mostly Staphylococcus aureus) that 
alters the bone remodelling process and that rapidly leads to severe bone loss, necrosis of the affected 
portion, and it may even spread to other parts of the body. Therefore it is one of the most appropriate 
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Figure 5: Optimally controlled solutions for the SIR model with therapies as controlled discrete inputs 
simulated for 15 days. Parameters are: S(0) = 0.3, 1(0) = 0.7, R(0) = 0,b = H = 0.02, j8 = 1800, V = 
100, p = 0.5, k = 50. Results have been computed by using a global optimal nonlinear solver. 



examples of multiscale infection, since it affects multiple biological scales. Indeed infection by S. au- 
reus starts at the intracellular level, but it subsequently involves also the tissue and the organ level. Bone 
remodelling (BR) is the cellular-level process by which the bone is continuously renewed as the result of 
an alternation of bone resorption, conducted by cells called osteoclasts, and bone formation, conducted 
by osteoblasts. They form together the so-called Basic Multicellular Unit (BMU), i.e. an ensemble of 
osteoclasts and osteoblasts that dissolve an area of the bone surface and then fill it. The BMU could be 
considered as an emerging scale intermediate between the cell and tissue levels. Osteomyelitis induces a 
severe inflammatory response followed by progressive bone destruction and loss of the vasculature and 
with a persistent chronic infection; this is further complicated by the rapid emergence of resistant strains 
of S. aureus. It has been shown that the infection prevents proliferation, induces apoptosis and inhibits 
mineralisation of cultured osteoblasts. Although effective treatment of this disease is very difficult, one 
of most used drug is the fusidic acid that acts as a bacteriostatic agent, and is usually combined with other 
antibiotics. Based on a recent work of the authors about the modelling of osteomyelitis and the compar- 
ison of different treatments |fT5l , we can define a hybrid version of that model and a treatment strategy 
combining antibiotic and anti-inflammatory therapies. The model describes how the bone remodelling 
dynamics of osteoblasts' (Ob) and osteoclasts' (Oc) population is affected by the S.aureus (B). 



O c 
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/t£12/(1+/i2 7)^22-/22 7 
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where a, and j3, are growth and death rates; gij describes the effectiveness of autocrine and paracrine 
regulation, i.e. the chemical interactions between osteoblasts and osteoclasts; fu models the impact of 
the infection on the autocrine and paracrine regulation; T\ is the discrete input modelling the antibiotic 
dosage; T2 models the dosage of a anti-inflammatory therapy; kj is the effectiveness of therapy T2; the 
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bacterial population follows a Gompertz curve with carrying capacity s and growth rate y#; and the output 
function y represents the bone density calculated as a function of Oc, Ob and their resorption (k\) and 
formation (#2) rates. Bone Mineral Density (BMD) measurement are often taken in medical practice, 
thus y is the quantity observed by the controller (a doctor who administers the treatment) after each visit. 

The optimal scheduling of bacteriostatic and anti-inflammatory therapies could be formulated with a 
control law similar to the SIR example. Since the combination of antibiotic and anti-inflammatory drugs 
has never been applied in medical practice, we could set the weight matrix of input therapies R so that 
Rn and R21 are much larger that Ru and 7?22- In this way, we give a high penalty to the combination of 
therapies Tl and T2 w.r.t. Tl and T2 taken individually. 



4 Discussion and conclusion 

The field of predictive models in biomedicine is challenged by the need of a better comprehension of 
the phenomenon of scales in the biological organisation, particularly their role in the transition between 
health and disease conditions. The multiscale modelling of molecules-cell-tissue-organ-body interac- 
tions is a key step in the process of identifying the most important parameters acting in a disease state 
and their calibration, linking basic research and clinical practice (therapies). Here we discuss how a de- 
scription based on hybrid dynamical systems is beneficial to the formulation of a multiscale modelling in 
biomedical processes. A second instance of the utility of hybrid dynamics approach stems from the pres- 
ence of multiple controls in biological systems. These controls could be framed as occurring naturally 
(the immune system response) or induced after the administration of a therapy. It is very often the case 
of multiple treatments, i.e. switching between therapies or combination of therapies which is the case 
considered in this study. In particular, we show the suitability of the D-CGF high-level process algebra 
to accomplish the task: here, its semantics is given in terms of a hybrid dynamical system where the 
large number of cell populations provides the basis for a continuous modelling approach and the dosage 
of multiple therapies is implemented as discrete switches. 

We believe that this approach opens several interesting directions in basic and clinical research. 
First, it is general and could be used for different tissues and organs or multiorgan diseases. Second, it is 
possible to extend the mathematical formulation to all the scales of biological organisation involved in 
the infection and recovery conditions. Third, the therapy could be complex and built on several nested 
and hierarchical protocols. One final comment: given the richness of examples provided by biological 
processes and medical therapies, there is basis for deriving interesting theories. We remind the bon mot 
of Stan Ulam: Ask not what mathematics can do for biology. Ask what biology can do for Mathematics. 
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